/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2007, 2008, 2009, 2010 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include "utilities.hpp"

#include <ql/time/schedule.hpp>
#include <ql/time/calendars/target.hpp>
#include <ql/quotes/simplequote.hpp>
#include <ql/instruments/europeanoption.hpp>
#include <ql/instruments/impliedvolatility.hpp>
#include <ql/processes/blackscholesprocess.hpp>
#include <ql/processes/hybridhestonhullwhiteprocess.hpp>
#include <ql/math/randomnumbers/rngtraits.hpp>
#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
#include <ql/math/optimization/simplex.hpp>
#include <ql/math/optimization/levenbergmarquardt.hpp>
#include <ql/math/statistics/generalstatistics.hpp>
#include <ql/math/statistics/sequencestatistics.hpp>
#include <ql/math/statistics/incrementalstatistics.hpp>
#include <ql/math/matrixutilities/svd.hpp>
#include <ql/time/daycounters/actual360.hpp>
#include <ql/time/daycounters/actual365fixed.hpp>
#include <ql/methods/montecarlo/multipathgenerator.hpp>
#include <ql/termstructures/yield/zerocurve.hpp>
#include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
#include <ql/models/equity/hestonmodelhelper.hpp>
#include <ql/models/shortrate/onefactormodels/hullwhite.hpp>
#include <ql/pricingengines/vanilla/analytichestonengine.hpp>
#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
#include <ql/pricingengines/vanilla/analytich1hwengine.hpp>
#include <ql/pricingengines/vanilla/mchestonhullwhiteengine.hpp>
#include <ql/pricingengines/vanilla/analyticbsmhullwhiteengine.hpp>
#include <ql/pricingengines/vanilla/analytichestonhullwhiteengine.hpp>
#include <ql/pricingengines/vanilla/fdhestonvanillaengine.hpp>
#include <ql/pricingengines/vanilla/fdhestonhullwhitevanillaengine.hpp>
#include <functional>

using namespace QuantLib;


TEST_CASE("HybridHestonHullWhiteProcess_BsmHullWhiteEngine", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing European option pricing for a BSM process"
                 " with one-factor Hull-White model...");

    SavedSettings backup;

    DayCounter dc = Actual365Fixed();

    const Date today = Date::todaysDate();
    const Date maturity = today + Period(20, Years);

    Settings::instance().evaluationDate() = today;

    const Handle<Quote> spot(
            std::make_shared<SimpleQuote>(100.0));
    std::shared_ptr < SimpleQuote > qRate = std::make_shared<SimpleQuote>(0.04);
    const Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
    std::shared_ptr < SimpleQuote > rRate = std::make_shared<SimpleQuote>(0.0525);
    const Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));
    std::shared_ptr < SimpleQuote > vol = std::make_shared<SimpleQuote>(0.25);
    const Handle<BlackVolTermStructure> volTS(flatVol(today, vol, dc));

    // FLOATING_POINT_EXCEPTION
    std::shared_ptr < HullWhite > hullWhiteModel =
            std::make_shared<HullWhite>(Handle<YieldTermStructure>(rTS), 0.00883, 0.00526);

    std::shared_ptr < BlackScholesMertonProcess > stochProcess =
            std::make_shared<BlackScholesMertonProcess>(spot, qTS, rTS, volTS);

    std::shared_ptr < Exercise > exercise = std::make_shared<EuropeanExercise>(maturity);

    Real fwd = spot->value() * qTS->discount(maturity) / rTS->discount(maturity);
    std::shared_ptr < StrikedTypePayoff > payoff =
            std::make_shared<PlainVanillaPayoff>(Option::Call, fwd);

    EuropeanOption option(payoff, exercise);

    const Real tol = 1e-8;
    const Real corr[] = {-0.75, -0.25, 0.0, 0.25, 0.75};
    const Volatility expectedVol[] = {0.217064577, 0.243995801,
                                      0.256402830, 0.268236596, 0.290461343};

    for (Size i = 0; i < LENGTH(corr); ++i) {
        std::shared_ptr < PricingEngine > bsmhwEngine =
                std::make_shared<AnalyticBSMHullWhiteEngine>(corr[i], stochProcess,
                                                             hullWhiteModel);

        option.setPricingEngine(bsmhwEngine);
        const Real npv = option.NPV();

        const Handle<BlackVolTermStructure> compVolTS(
                flatVol(today, expectedVol[i], dc));

        std::shared_ptr < BlackScholesMertonProcess > bsProcess =
                std::make_shared<BlackScholesMertonProcess>(spot, qTS, rTS, compVolTS);
        std::shared_ptr < PricingEngine > bsEngine =
                std::make_shared<AnalyticEuropeanEngine>(bsProcess);

        EuropeanOption comp(payoff, exercise);
        comp.setPricingEngine(bsEngine);

        Volatility impliedVol =
                comp.impliedVolatility(npv, bsProcess, 1e-10, 100);

        if (std::fabs(impliedVol - expectedVol[i]) > tol) {
            FAIL("Failed to reproduce implied volatility"
                         << "\n    calculated: " << impliedVol
                         << "\n    expected  : " << expectedVol[i]);
        }
        if (std::fabs((comp.NPV() - npv) / npv) > tol) {
            FAIL("Failed to reproduce NPV"
                         << "\n    calculated: " << npv
                         << "\n    expected  : " << comp.NPV());
        }
        if (std::fabs(comp.delta() - option.delta()) > tol) {
            FAIL("Failed to reproduce NPV"
                         << "\n    calculated: " << npv
                         << "\n    expected  : " << comp.NPV());
        }
        if (std::fabs((comp.gamma() - option.gamma()) / npv) > tol) {
            FAIL("Failed to reproduce NPV"
                         << "\n    calculated: " << npv
                         << "\n    expected  : " << comp.NPV());
        }
        if (std::fabs((comp.theta() - option.theta()) / npv) > tol) {
            FAIL("Failed to reproduce NPV"
                         << "\n    calculated: " << npv
                         << "\n    expected  : " << comp.NPV());
        }
        if (std::fabs((comp.vega() - option.vega()) / npv) > tol) {
            FAIL("Failed to reproduce NPV"
                         << "\n    calculated: " << npv
                         << "\n    expected  : " << comp.NPV());
        }
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_CompareBsmHWandHestonHW", "[HybridHestonHullWhiteProcess]") {
    INFO("Comparing European option pricing for a BSM process"
                 " with one-factor Hull-White model...");

    SavedSettings backup;

    DayCounter dc = Actual365Fixed();

    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    const Handle<Quote> spot(
            std::make_shared<SimpleQuote>(100.0));
    std::vector<Date> dates;
    std::vector<Rate> rates, divRates;

    for (Size i = 0; i <= 40; ++i) {
        dates.emplace_back(today + Period(i, Years));
        // FLOATING_POINT_EXCEPTION
        rates.emplace_back(0.01 + 0.0002 * std::exp(std::sin(i / 4.0)));
        divRates.emplace_back(0.02 + 0.0001 * std::exp(std::sin(i / 5.0)));
    }

    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100));
    const Handle<YieldTermStructure> rTS(
            std::make_shared<ZeroCurve>(dates, rates, dc));
    const Handle<YieldTermStructure> qTS(
            std::make_shared<ZeroCurve>(dates, divRates, dc));

    std::shared_ptr < SimpleQuote > vol = std::make_shared<SimpleQuote>(0.25);
    const Handle<BlackVolTermStructure> volTS(flatVol(today, vol, dc));

    std::shared_ptr < BlackScholesMertonProcess > bsmProcess =
            std::make_shared<BlackScholesMertonProcess>(spot, qTS, rTS, volTS);

    std::shared_ptr < HestonProcess > hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, spot,
                                            vol->value() * vol->value(), 1.0,
                                            vol->value() * vol->value(), 1e-4, 0.0);

    std::shared_ptr < HestonModel > hestonModel = std::make_shared<HestonModel>(hestonProcess);

    std::shared_ptr < HullWhite > hullWhiteModel =
            std::make_shared<HullWhite>(Handle<YieldTermStructure>(rTS), 0.01, 0.01);

    std::shared_ptr < PricingEngine > bsmhwEngine =
            std::make_shared<AnalyticBSMHullWhiteEngine>(0.0, bsmProcess, hullWhiteModel);

    std::shared_ptr < PricingEngine > hestonHwEngine =
            std::make_shared<AnalyticHestonHullWhiteEngine>(hestonModel, hullWhiteModel, 128);


    const Real tol = 1e-5;
    const Real strike[] = {0.25, 0.5, 0.75, 0.8, 0.9,
                           1.0, 1.1, 1.2, 1.5, 2.0, 4.0};
    const Size maturity[] = {1, 2, 3, 5, 10, 15, 20, 25, 30};
    const Option::Type types[] = {Option::Put, Option::Call};

    for (Size i = 0; i < LENGTH(types); ++i) {
        for (Size j = 0; j < LENGTH(strike); ++j) {
            for (Size l = 0; l < LENGTH(maturity); ++l) {
                const Date maturityDate = today + Period(maturity[l], Years);

                std::shared_ptr < Exercise > exercise =
                        std::make_shared<EuropeanExercise>(maturityDate);

                Real fwd = strike[j] * spot->value()
                           * qTS->discount(maturityDate) / rTS->discount(maturityDate);

                std::shared_ptr < StrikedTypePayoff > payoff = std::make_shared<PlainVanillaPayoff>(types[i], fwd);

                EuropeanOption option(payoff, exercise);

                option.setPricingEngine(bsmhwEngine);
                const Real calculated = option.NPV();

                option.setPricingEngine(hestonHwEngine);
                const Real expected = option.NPV();

                if (std::fabs(calculated - expected) > calculated * tol &&
                    std::fabs(calculated - expected) > tol) {
                    FAIL_CHECK("Failed to reproduce npvs"
                                       << "\n    calculated: " << calculated
                                       << "\n    expected  : " << expected
                                       << "\n    strike    : " << strike[j]
                                       << "\n    maturity  : " << maturity[l]
                                       << "\n    type      : "
                                       << ((types[i] == Option::Put) ? "Put" : "Call")
                    );
                }
            }
        }
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_ZeroBondPricing", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing Monte-Carlo zero bond pricing...");

    SavedSettings backup;

    DayCounter dc = Actual360();
    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    // construct a strange yield curve to check drifts and discounting
    // of the joint stochastic process

    std::vector<Date> dates;
    std::vector<Time> times;
    std::vector<Rate> rates;

    dates.emplace_back(today);
    rates.emplace_back(0.02);
    times.emplace_back(0.0);
    for (Size i = 120; i < 240; ++i) {
        dates.emplace_back(today + Period(i, Months));
        rates.emplace_back(0.02 + 0.0002 * std::exp(std::sin(i / 8.0)));
        times.emplace_back(dc.yearFraction(today, dates.back()));
    }

    const Date maturity = dates.back() + Period(10, Years);
    dates.emplace_back(maturity);
    rates.emplace_back(0.04);
    times.emplace_back(dc.yearFraction(today, dates.back()));

    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100));

    const Handle<YieldTermStructure> ts(
            std::make_shared<ZeroCurve>(dates, rates, dc));
    const Handle<YieldTermStructure> ds(flatRate(today, 0.0, dc));

    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(ts, ds, s0, 0.02, 1.0, 0.2, 0.5, -0.8);
    const std::shared_ptr<HullWhiteForwardProcess> hwProcess =
            std::make_shared<HullWhiteForwardProcess>(ts, 0.05, 0.05);
    hwProcess->setForwardMeasureTime(dc.yearFraction(today, maturity));
    const std::shared_ptr<HullWhite> hwModel = std::make_shared<HullWhite>(ts, 0.05, 0.05);

    const std::shared_ptr<HybridHestonHullWhiteProcess> jointProcess =
            std::make_shared<HybridHestonHullWhiteProcess>(hestonProcess, hwProcess, -0.4);

    TimeGrid grid(times.begin(), times.end() - 1);

    typedef SobolBrownianBridgeRsg rsg_type;
    typedef MultiPathGenerator<rsg_type>::sample_type sample_type;

    const Size factors = jointProcess->factors();
    const Size steps = grid.size() - 1;
    rsg_type rsg = rsg_type(factors, steps);
    MultiPathGenerator<rsg_type> generator(jointProcess, grid, rsg, false);

    const Size m = 90;
    std::vector<GeneralStatistics> zeroStat(m);
    std::vector<GeneralStatistics> optionStat(m);

    const Size nrTrails = 8191;
    const Size optionTenor = 24;
    const DiscountFactor strike = 0.5;

    for (Size i = 0; i < nrTrails; ++i) {
        sample_type path = generator.next();

        for (Size j = 1; j < m; ++j) {
            const Time t = grid[j];            // zero end and option maturity
            const Time T = grid[j + optionTenor];// maturity of zero bond
            // of option

            Array states(3);
            Array optionStates(3);
            for (Size k = 0; k < jointProcess->size(); ++k) {
                states[k] = path.value[k][j];
                optionStates[k] = path.value[k][j + optionTenor];
            }

            const DiscountFactor zeroBond
                    = 1.0 / jointProcess->numeraire(t, states);
            const DiscountFactor zeroOption = zeroBond
                                              * std::max(0.0, hwModel->discountBond(t, T, states[2]) - strike);

            zeroStat[j].add(zeroBond);
            optionStat[j].add(zeroOption);
        }
    }

    for (Size j = 1; j < m; ++j) {
        const Time t = grid[j];
        Real calculated = zeroStat[j].mean();
        Real expected = ts->discount(t);

        if (std::fabs(calculated - expected) > 0.03) {
            FAIL_CHECK("Failed to reproduce expected zero bond prices"
                               << "\n   t:          " << t
                               << "\n   calculated: " << calculated
                               << "\n   expected:   " << expected);
        }

        const Time T = grid[j + optionTenor];

        calculated = optionStat[j].mean();
        expected = hwModel->discountBondOption(Option::Call, strike, t, T);

        if (std::fabs(calculated - expected) > 0.0035) {
            FAIL_CHECK("Failed to reproduce expected zero bond option prices"
                               << "\n   t:          " << t
                               << "\n   T:          " << T
                               << "\n   calculated: " << calculated
                               << "\n   expected:   " << expected);
        }
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_McVanillaPricing", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing Monte-Carlo vanilla option pricing...");

    SavedSettings backup;

    DayCounter dc = Actual360();
    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    // construct a strange yield curve to check drifts and discounting
    // of the joint stochastic process

    std::vector<Date> dates;
    std::vector<Rate> rates, divRates;

    for (Size i = 0; i <= 40; ++i) {
        dates.emplace_back(today + Period(i, Years));
        // FLOATING_POINT_EXCEPTION
        rates.emplace_back(0.03 + 0.0003 * std::exp(std::sin(i / 4.0)));
        divRates.emplace_back(0.02 + 0.0001 * std::exp(std::sin(i / 5.0)));
    }

    const Date maturity = today + Period(20, Years);

    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100));
    const Handle<YieldTermStructure> rTS(
            std::make_shared<ZeroCurve>(dates, rates, dc));
    const Handle<YieldTermStructure> qTS(
            std::make_shared<ZeroCurve>(dates, divRates, dc));
    std::shared_ptr < SimpleQuote > vol = std::make_shared<SimpleQuote>(0.25);
    const Handle<BlackVolTermStructure> volTS(flatVol(today, vol, dc));

    const std::shared_ptr<BlackScholesMertonProcess> bsmProcess =
            std::make_shared<BlackScholesMertonProcess>(s0, qTS, rTS, volTS);
    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, s0, 0.0625, 0.5, 0.0625, 1e-5, 0.3);
    const std::shared_ptr<HullWhiteForwardProcess> hwProcess =
            std::make_shared<HullWhiteForwardProcess>(rTS, 0.01, 0.01);
    hwProcess->setForwardMeasureTime(dc.yearFraction(today, maturity));

    const Real tol = 0.05;
    const Real corr[] = {-0.9, -0.5, 0.0, 0.5, 0.9};
    const Real strike[] = {100};

    for (Size i = 0; i < LENGTH(corr); ++i) {
        for (Size j = 0; j < LENGTH(strike); ++j) {
            std::shared_ptr < HybridHestonHullWhiteProcess > jointProcess =
                    std::make_shared<HybridHestonHullWhiteProcess>(hestonProcess,
                                                                   hwProcess, corr[i]);

            std::shared_ptr < StrikedTypePayoff > payoff =
                    std::make_shared<PlainVanillaPayoff>(Option::Put, strike[j]);
            std::shared_ptr < Exercise > exercise =
                    std::make_shared<EuropeanExercise>(maturity);

            VanillaOption optionHestonHW(payoff, exercise);
            std::shared_ptr < PricingEngine > engine =
                    MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
                            .withSteps(1)
                            .withAntitheticVariate()
                            .withControlVariate()
                            .withAbsoluteTolerance(tol)
                            .withSeed(42);

            optionHestonHW.setPricingEngine(engine);

            const std::shared_ptr<HullWhite> hwModel =
                    std::make_shared<HullWhite>(Handle<YieldTermStructure>(rTS),
                                                hwProcess->a(), hwProcess->sigma());

            VanillaOption optionBsmHW(payoff, exercise);
            optionBsmHW.setPricingEngine(
                    std::make_shared<AnalyticBSMHullWhiteEngine>(corr[i], bsmProcess,
                                                                 hwModel));

            const Real calculated = optionHestonHW.NPV();
            const Real error = optionHestonHW.errorEstimate();
            const Real expected = optionBsmHW.NPV();

            if ((corr[i] != 0.0 && std::fabs(calculated - expected) > 3 * error)
                || (corr[i] == 0.0 && std::fabs(calculated - expected) > 1e-4)) {
                FAIL_CHECK("Failed to reproduce BSM-HW vanilla prices"
                                   << "\n   corr:       " << corr[i]
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   error:      " << error
                                   << "\n   expected:   " << expected);
            }
        }
    }
}


TEST_CASE("HybridHestonHullWhiteProcess_McPureHestonPricing", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing Monte-Carlo Heston option pricing...");

    SavedSettings backup;

    DayCounter dc = Actual360();
    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    // construct a strange yield curve to check drifts and discounting
    // of the joint stochastic process

    std::vector<Date> dates;
    std::vector<Rate> rates, divRates;

    for (Size i = 0; i <= 100; ++i) {
        dates.emplace_back(today + Period(i, Months));
        // FLOATING_POINT_EXCEPTION
        rates.emplace_back(0.02 + 0.0002 * std::exp(std::sin(i / 10.0)));
        divRates.emplace_back(0.02 + 0.0001 * std::exp(std::sin(i / 20.0)));
    }

    const Date maturity = today + Period(2, Years);

    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100));
    const Handle<YieldTermStructure> rTS(
            std::make_shared<ZeroCurve>(dates, rates, dc));
    const Handle<YieldTermStructure> qTS(
            std::make_shared<ZeroCurve>(dates, divRates, dc));

    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, s0, 0.08, 1.5, 0.0625, 0.5, -0.8);
    const std::shared_ptr<HullWhiteForwardProcess> hwProcess =
            std::make_shared<HullWhiteForwardProcess>(rTS, 0.1, 1e-8);
    hwProcess->setForwardMeasureTime(dc.yearFraction(
            today, maturity + Period(1, Years)));

    const Real tol = 0.001;
    const Real corr[] = {-0.45, 0.45, 0.25};
    const Real strike[] = {100, 75, 50, 150};

    for (Size i = 0; i < LENGTH(corr); ++i) {
        for (Size j = 0; j < LENGTH(strike); ++j) {
            std::shared_ptr < HybridHestonHullWhiteProcess > jointProcess =
                    std::make_shared<HybridHestonHullWhiteProcess>(
                            hestonProcess, hwProcess,
                            corr[i], HybridHestonHullWhiteProcess::Euler);

            std::shared_ptr < StrikedTypePayoff > payoff =
                    std::make_shared<PlainVanillaPayoff>(Option::Put, strike[j]);
            std::shared_ptr < Exercise > exercise =
                    std::make_shared<EuropeanExercise>(maturity);

            VanillaOption optionHestonHW(payoff, exercise);
            VanillaOption optionPureHeston(payoff, exercise);
            optionPureHeston.setPricingEngine(
                    std::make_shared<AnalyticHestonEngine>(
                            std::make_shared<HestonModel>(hestonProcess)));

            Real expected = optionPureHeston.NPV();

            optionHestonHW.setPricingEngine(
                    MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
                            .withSteps(2)
                            .withAntitheticVariate()
                            .withControlVariate()
                            .withAbsoluteTolerance(tol)
                            .withSeed(42));

            Real calculated = optionHestonHW.NPV();
            Real error = optionHestonHW.errorEstimate();

            if (std::fabs(calculated - expected) > 3 * error
                && std::fabs(calculated - expected) > tol) {
                FAIL_CHECK("Failed to reproduce pure heston vanilla prices"
                                   << "\n   corr:       " << corr[i]
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   error:      " << error
                                   << "\n   expected:   " << expected);
            }
        }
    }
}


TEST_CASE("HybridHestonHullWhiteProcess_AnalyticHestonHullWhitePricing", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing analytic Heston Hull-White option pricing...");

    SavedSettings backup;

    DayCounter dc = Actual360();
    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    // construct a strange yield curve to check drifts and discounting
    // of the joint stochastic process

    std::vector<Date> dates;
    std::vector<Rate> rates, divRates;

    for (Size i = 0; i <= 40; ++i) {
        dates.emplace_back(today + Period(i, Years));
        // FLOATING_POINT_EXCEPTION
        rates.emplace_back(0.03 + 0.0001 * std::exp(std::sin(i / 4.0)));
        divRates.emplace_back(0.02 + 0.0002 * std::exp(std::sin(i / 3.0)));
    }

    const Date maturity = today + Period(5, Years);
    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100));
    const Handle<YieldTermStructure> rTS(
            std::make_shared<ZeroCurve>(dates, rates, dc));
    const Handle<YieldTermStructure> qTS(
            std::make_shared<ZeroCurve>(dates, divRates, dc));

    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, s0, 0.08, 1.5, 0.0625, 0.5, -0.8);
    const std::shared_ptr<HestonModel> hestonModel =
            std::make_shared<HestonModel>(hestonProcess);

    const std::shared_ptr<HullWhiteForwardProcess> hwFwdProcess =
            std::make_shared<HullWhiteForwardProcess>(rTS, 0.01, 0.01);
    hwFwdProcess->setForwardMeasureTime(dc.yearFraction(today, maturity));
    const std::shared_ptr<HullWhite> hullWhiteModel = std::make_shared<HullWhite>(
            rTS, hwFwdProcess->a(), hwFwdProcess->sigma());

    const Real tol = 0.002;
    const Real strike[] = {80, 120};
    const Option::Type types[] = {Option::Put, Option::Call};

    for (Size i = 0; i < LENGTH(types); ++i) {
        for (Size j = 0; j < LENGTH(strike); ++j) {
            std::shared_ptr < HybridHestonHullWhiteProcess > jointProcess =
                    std::make_shared<HybridHestonHullWhiteProcess>(
                            hestonProcess, hwFwdProcess, 0.0,
                            HybridHestonHullWhiteProcess::Euler);

            std::shared_ptr < StrikedTypePayoff > payoff =
                    std::make_shared<PlainVanillaPayoff>(types[i], strike[j]);
            std::shared_ptr < Exercise > exercise =
                    std::make_shared<EuropeanExercise>(maturity);

            VanillaOption optionHestonHW(payoff, exercise);
            optionHestonHW.setPricingEngine(
                    MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
                            .withSteps(1)
                            .withAntitheticVariate()
                            .withControlVariate()
                            .withAbsoluteTolerance(tol)
                            .withSeed(42));

            VanillaOption optionPureHeston(payoff, exercise);
            optionPureHeston.setPricingEngine(
                    std::make_shared<AnalyticHestonHullWhiteEngine>(hestonModel,
                                                                    hullWhiteModel, 128));

            Real calculated = optionHestonHW.NPV();
            Real error = optionHestonHW.errorEstimate();
            Real expected = optionPureHeston.NPV();

            if (std::fabs(calculated - expected) > 3 * error
                && std::fabs(calculated - expected) > tol) {
                FAIL_CHECK("Failed to reproduce hw heston vanilla prices"
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   error:      " << error
                                   << "\n   expected:   " << expected);
            }
        }
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_CallableEquityPricing", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing the pricing of a callable equity product...");

    SavedSettings backup;

    /*
       For the definition of the example product see
       Alexander Giese, On the Pricing of Auto-Callable Equity
       Structures in the Presence of Stochastic Volatility and
       Stochastic Interest Rates .
       http://workshop.mathfinance.de/2006/papers/giese/slides.pdf
    */

    const Size maturity = 7;
    DayCounter dc = Actual365Fixed();
    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    Handle<Quote> spot(std::make_shared<SimpleQuote>(100.0));
    std::shared_ptr < SimpleQuote > qRate = std::make_shared<SimpleQuote>(0.04);
    Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
    std::shared_ptr < SimpleQuote > rRate = std::make_shared<SimpleQuote>(0.04);
    Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));

    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, spot, 0.0625, 1.0,
                                            0.24 * 0.24, 1e-4, 0.0);
    // FLOATING_POINT_EXCEPTION
    const std::shared_ptr<HullWhiteForwardProcess> hwProcess =
            std::make_shared<HullWhiteForwardProcess>(rTS, 0.00883, 0.00526);
    hwProcess->setForwardMeasureTime(
            dc.yearFraction(today, today + Period(maturity + 1, Years)));

    const std::shared_ptr<HybridHestonHullWhiteProcess> jointProcess =
            std::make_shared<HybridHestonHullWhiteProcess>(hestonProcess, hwProcess, -0.4);

    Schedule schedule(today, today + Period(maturity, Years),
                      Period(1, Years), TARGET(),
                      Following, Following,
                      DateGeneration::Forward, false);

    std::vector<Time> times(maturity + 1);
    std::transform(schedule.begin(), schedule.end(), times.begin(),
                   [dc, today](const Date &d) { return dc.yearFraction(today, d, Date(), Date()); });

    for (Size i = 0; i <= maturity; ++i)
        times[i] = static_cast<Time>(i);

    TimeGrid grid(times.begin(), times.end());

    std::vector<Real> redemption(maturity);
    for (Size i = 0; i < maturity; ++i) {
        redemption[i] = 1.07 + 0.03 * i;
    }

    typedef PseudoRandom::rsg_type rsg_type;
    typedef MultiPathGenerator<rsg_type>::sample_type sample_type;

    BigNatural seed = 42;
    rsg_type rsg = PseudoRandom::make_sequence_generator(
            jointProcess->factors() * (grid.size() - 1), seed);

    MultiPathGenerator<rsg_type> generator(jointProcess, grid, rsg, false);
    GeneralStatistics stat;

    Real antitheticPayoff = 0;
    const Size nrTrails = 40000;
    for (Size i = 0; i < nrTrails; ++i) {
        const bool antithetic = (i % 2) == 0 ? false : true;

        sample_type path = antithetic ? generator.antithetic()
                                      : generator.next();

        Real payoff = 0;
        for (Size j = 1; j <= maturity; ++j) {
            if (path.value[0][j] > spot->value()) {
                Array states(3);
                for (Size k = 0; k < 3; ++k) {
                    states[k] = path.value[k][j];
                }
                payoff = redemption[j - 1]
                         / jointProcess->numeraire(grid[j], states);
                break;
            } else if (j == maturity) {
                Array states(3);
                for (Size k = 0; k < 3; ++k) {
                    states[k] = path.value[k][j];
                }
                payoff = 1.0 / jointProcess->numeraire(grid[j], states);
            }
        }

        if (antithetic) {
            stat.add(0.5 * (antitheticPayoff + payoff));
        } else {
            antitheticPayoff = payoff;
        }
    }

    const Real expected = 0.938;
    const Real calculated = stat.mean();
    const Real error = stat.errorEstimate();

    if (std::fabs(expected - calculated) > 3 * error) {
        FAIL_CHECK("Failed to reproduce auto-callable equity structure price"
                           << "\n   calculated: " << calculated
                           << "\n   error:      " << error
                           << "\n   expected:   " << expected);
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_DiscretizationError", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing the discretization error of the "
                 "Heston Hull-White process...");

    SavedSettings backup;

    DayCounter dc = Actual360();
    const Date today = Date::todaysDate();

    Settings::instance().evaluationDate() = today;

    // construct a strange yield curve to check drifts and discounting
    // of the joint stochastic process

    std::vector<Date> dates;
    std::vector<Rate> rates, divRates;

    for (Size i = 0; i <= 31; ++i) {
        dates.emplace_back(today + Period(i, Years));
        // FLOATING_POINT_EXCEPTION
        rates.emplace_back(0.04 + 0.0001 * std::exp(std::sin(double(i))));
        divRates.emplace_back(0.04 + 0.0001 * std::exp(std::sin(double(i))));
    }

    const Date maturity = today + Period(10, Years);
    const Volatility v = 0.25;

    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100));
    const std::shared_ptr<SimpleQuote> vol = std::make_shared<SimpleQuote>(v);
    const Handle<BlackVolTermStructure> volTS(flatVol(today, vol, dc));
    const Handle<YieldTermStructure> rTS(
            std::make_shared<ZeroCurve>(dates, rates, dc));
    const Handle<YieldTermStructure> qTS(
            std::make_shared<ZeroCurve>(dates, divRates, dc));

    const std::shared_ptr<BlackScholesMertonProcess> bsmProcess =
            std::make_shared<BlackScholesMertonProcess>(s0, qTS, rTS, volTS);

    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, s0, v * v, 1, v * v, 1e-6, -0.4);

    const std::shared_ptr<HullWhiteForwardProcess> hwProcess =
            std::make_shared<HullWhiteForwardProcess>(rTS, 0.01, 0.01);
    hwProcess->setForwardMeasureTime(20.1472222222222222);

    const Real tol = 0.05;
    const Real corr[] = {-0.85, 0.5};
    const Real strike[] = {50, 100, 125};

    for (Size i = 0; i < LENGTH(corr); ++i) {
        for (Size j = 0; j < LENGTH(strike); ++j) {
            std::shared_ptr < StrikedTypePayoff > payoff =
                    std::make_shared<PlainVanillaPayoff>(Option::Put, strike[j]);
            std::shared_ptr < Exercise > exercise =
                    std::make_shared<EuropeanExercise>(maturity);

            VanillaOption optionBsmHW(payoff, exercise);
            const std::shared_ptr<HullWhite> hwModel = std::make_shared<HullWhite>(
                    rTS, hwProcess->a(), hwProcess->sigma());
            optionBsmHW.setPricingEngine(
                    std::make_shared<AnalyticBSMHullWhiteEngine>(corr[i], bsmProcess,
                                                                 hwModel));

            Real expected = optionBsmHW.NPV();

            VanillaOption optionHestonHW(payoff, exercise);
            std::shared_ptr < HybridHestonHullWhiteProcess > jointProcess =
                    std::make_shared<HybridHestonHullWhiteProcess>(hestonProcess,
                                                                   hwProcess, corr[i]);
            optionHestonHW.setPricingEngine(
                    MakeMCHestonHullWhiteEngine<PseudoRandom>(jointProcess)
                            .withSteps(1)
                            .withAntitheticVariate()
                            .withAbsoluteTolerance(tol)
                            .withSeed(42));

            Real calculated = optionHestonHW.NPV();
            Real error = optionHestonHW.errorEstimate();

            if ((std::fabs(calculated - expected) > 3 * error
                 && std::fabs(calculated - expected) > 1e-5)) {
                FAIL_CHECK("Failed to reproduce discretization error"
                                   << "\n   corr:       " << corr[i]
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   error:      " << error
                                   << "\n   expected:   " << expected);
            }
        }
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_FdmHestonHullWhiteEngine", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing the FDM Heston Hull-White engine...");

    SavedSettings backup;

    const Date today = Date(28, March, 2004);
    Settings::instance().evaluationDate() = today;
    const Date exerciseDate = Date(28, March, 2012);
    DayCounter dc = Actual365Fixed();

    Handle<Quote> s0(std::make_shared<SimpleQuote>(100.0));

    const Handle<YieldTermStructure> rTS(flatRate(0.05, dc));
    const Handle<YieldTermStructure> qTS(flatRate(0.02, dc));

    const Volatility vol = 0.30;
    const Handle<BlackVolTermStructure> volTS(flatVol(vol, dc));

    const Real v0 = vol * vol;
    std::shared_ptr < HestonProcess > hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, s0, v0, 1.0, v0, 0.000001, 0.0);

    std::shared_ptr < BlackScholesMertonProcess > stochProcess =
            std::make_shared<BlackScholesMertonProcess>(s0, qTS, rTS, volTS);

    std::shared_ptr < HullWhiteProcess > hwProcess =
            std::make_shared<HullWhiteProcess>(rTS, 0.00883, 0.01);
    std::shared_ptr < HullWhite > hwModel =
            std::make_shared<HullWhite>(rTS, hwProcess->a(), hwProcess->sigma());

    std::shared_ptr < Exercise > exercise = std::make_shared<EuropeanExercise>(exerciseDate);
    const Real corr[] = {-0.85, 0.5};
    const Real strike[] = {75, 120, 160};

    for (Size i = 0; i < LENGTH(corr); ++i) {
        for (Size j = 0; j < LENGTH(strike); ++j) {
            std::shared_ptr < StrikedTypePayoff > payoff =
                    std::make_shared<PlainVanillaPayoff>(Option::Call, strike[j]);
            VanillaOption option(payoff, exercise);

            option.setPricingEngine(std::make_shared<FdHestonHullWhiteVanillaEngine>(
                    std::make_shared<HestonModel>(hestonProcess),
                    hwProcess, corr[i], 50, 200, 10, 15));
            const Real calculated = option.NPV();
            const Real calculatedDelta = option.delta();
            const Real calculatedGamma = option.gamma();

            option.setPricingEngine(std::make_shared<AnalyticBSMHullWhiteEngine>(corr[i],
                                                                                 stochProcess, hwModel));
            const Real expected = option.NPV();
            const Real expectedDelta = option.delta();
            const Real expectedGamma = option.gamma();

            const Real npvTol = 0.01;
            if (std::fabs(calculated - expected) > npvTol) {
                FAIL_CHECK("Failed to reproduce analytic npv values"
                                   << "\n   corr:       " << corr[i]
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   expected:   " << expected);
            }
            const Real deltaTol = 0.001;
            if (std::fabs(calculatedDelta - expectedDelta) > deltaTol) {
                FAIL_CHECK("Failed to reproduce analytic delta values"
                                   << "\n   corr:       " << corr[i]
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   expected:   " << expected);
            }
            const Real gammaTol = 0.001;
            if (std::fabs(calculatedGamma - expectedGamma) > gammaTol) {
                FAIL_CHECK("Failed to reproduce analytic gamma values"
                                   << "\n   corr:       " << corr[i]
                                   << "\n   strike:     " << strike[j]
                                   << "\n   calculated: " << calculated
                                   << "\n   expected:   " << expected);
            }
        }
    }
}


namespace {

    struct HestonModelData {
        const char *const name;
        Real v0;
        Real kappa;
        Real theta;
        Real sigma;
        Real rho;
        Real r;
        Real q;
    };

    HestonModelData hestonModels[] = {
            // ADI finite difference schemes for option pricing in the
            // Heston model with correlation, K.J. in t'Hout and S. Foulon,
            {"'t Hout case 1",   0.04,   1.5,    0.04,   0.3,    -0.9,    0.025,  0.0},
            {"'t Hout case 2",   0.12,   3.0,    0.12,   0.04,   0.6,     0.01,   0.04},
            {"'t Hout case 3",   0.0707, 0.6067, 0.0707, 0.2928, -0.7571, 0.03,   0.0},
            {"'t Hout case 4",   0.06,   2.5,    0.06,   0.5,    -0.1,    0.0507, 0.0469},
            // Efficient numerical methods for pricing American options under
            // stochastic volatility, Samuli Ikonen and Jari Toivanen,
            {"Ikonen-Toivanen",  0.0625, 5,      0.16,   0.9,    0.1,     0.1,    0.0},
            // Not-so-complex logarithms in the Heston model,
            // Christian Kahl and Peter JГ¤ckel
            {"Kahl-Jaeckel",     0.16,   1.0,    0.16,   2.0,    -0.8,    0.0,    0.0},
            // self defined test cases
            {"Equity case",      0.07,   2.0,    0.04,   0.55,   -0.8,    0.03,   0.035},
            {"high correlation", 0.07,   1.0,    0.04,   0.55,   0.995,   0.02,   0.04},
            {"low Vol-Of-Vol",   0.07,   1.0,    0.04,   0.001,  -0.75,   0.04,   0.03},
            {"kappaEqSigRho",    0.07,   0.4,    0.04,   0.5,    0.8,     0.03,   0.03}
    };

    struct HullWhiteModelData {
        const char *const name;
        Real a;
        Real sigma;
    };

    HullWhiteModelData hullWhiteModels[] = {
            {"EUR-2003", 0.00883, 0.00631}
    };


    struct SchemeData {
        const char *const name;
        FdmSchemeDesc schemeDesc;
    };

    SchemeData schemes[] = {
            {"HV2", FdmSchemeDesc::Hundsdorfer()},
            {"HV1", FdmSchemeDesc::ModifiedHundsdorfer()},
            {"CS",  FdmSchemeDesc::CraigSneyd()},
            {"MCS", FdmSchemeDesc::ModifiedCraigSneyd()},
            {"DS",  FdmSchemeDesc::Douglas()}
    };

    struct VanillaOptionData {
        Real strike;
        Time maturity;
        Option::Type optionType;
    };

    std::shared_ptr<HestonProcess> makeHestonProcess(
            const HestonModelData &params) {

        Handle<Quote> spot(
                std::make_shared<SimpleQuote>(100));

        DayCounter dayCounter = Actual365Fixed();
        Handle<YieldTermStructure> rTS(flatRate(params.r, dayCounter));
        Handle<YieldTermStructure> qTS(flatRate(params.q, dayCounter));

        return std::make_shared<HestonProcess>(rTS, qTS, spot, params.v0, params.kappa, params.theta, params.sigma,
                                               params.rho);
    }

    std::shared_ptr<VanillaOption> makeVanillaOption(
            const VanillaOptionData &params) {

        Date maturity = Date(Settings::instance().evaluationDate())
                        + Period(Size(params.maturity * 365), Days);
        std::shared_ptr < Exercise > exercise = std::make_shared<EuropeanExercise>(maturity);
        std::shared_ptr < StrikedTypePayoff > payoff =
                std::make_shared<PlainVanillaPayoff>(params.optionType, params.strike);

        return std::make_shared<VanillaOption>(payoff, exercise);
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_BsmHullWhitePricing", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing convergence speed of Heston-Hull-White engine...");

    SavedSettings backup;

    Date today(27, December, 2004);
    Settings::instance().evaluationDate() = today;

    Real maturity = 5.0;
    Real equityIrCorr = -0.4;
    Real strikes[] = {75, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 140, 150};
    Size listOfTimeStepsPerYear[] = {20};

    HestonModelData hestonModelData
            = {"BSM-HW Model", 0.09, 1.0, 0.09, QL_EPSILON, 0.0, 0.04, 0.03};
    HullWhiteModelData hwModelData = hullWhiteModels[0];
    bool controlVariate[] = {true, false};

    std::shared_ptr < HestonProcess > hp(makeHestonProcess(hestonModelData));
    std::shared_ptr < HestonModel > hestonModel = std::make_shared<HestonModel>(hp);

    std::shared_ptr < HullWhiteProcess > hwProcess =
            std::make_shared<HullWhiteProcess>(hp->riskFreeRate(),
                                               hwModelData.a, hwModelData.sigma);
    std::shared_ptr < HullWhite > hullWhiteModel =
            std::make_shared<HullWhite>(hp->riskFreeRate(),
                                        hwProcess->a(), hwProcess->sigma());


    std::shared_ptr < BlackScholesMertonProcess > bsmProcess =
            std::make_shared<BlackScholesMertonProcess>(
                    hp->s0(), hp->dividendYield(), hp->riskFreeRate(),
                    Handle<BlackVolTermStructure>(
                            flatVol(today, std::sqrt(hestonModelData.theta),
                                    hp->riskFreeRate()->dayCounter())));

    std::shared_ptr < PricingEngine > bsmhwEngine =
            std::make_shared<AnalyticBSMHullWhiteEngine>(equityIrCorr, bsmProcess,
                                                         hullWhiteModel);

    Real tolWithCV[] = {2e-4, 2e-4, 2e-4, 2e-4, 0.01};
    Real tolWithOutCV[] = {5e-3, 5e-3, 5e-3, 5e-3, 0.02};
    for (Size l = 0; l < LENGTH(schemes); ++l) {
        SchemeData scheme = schemes[l];
        for (Size i = 0; i < LENGTH(controlVariate); ++i) {
            for (Size u = 0; u < LENGTH(listOfTimeStepsPerYear); ++u) {
                Size tSteps = Size(maturity * listOfTimeStepsPerYear[u]);

                std::shared_ptr < FdHestonHullWhiteVanillaEngine > fdEngine =
                        std::make_shared<FdHestonHullWhiteVanillaEngine>(
                                hestonModel, hwProcess, equityIrCorr,
                                tSteps, 400, 2, 10, 0, controlVariate[i],
                                scheme.schemeDesc);
                fdEngine->enableMultipleStrikesCaching(
                        std::vector<Real>(strikes, strikes + LENGTH(strikes)));

                Real avgPriceDiff = 0.0;
                for (Size k = 0; k < LENGTH(strikes); ++k) {
                    VanillaOptionData optionData
                            = {strikes[k], maturity, Option::Call};
                    std::shared_ptr < VanillaOption > option
                            = makeVanillaOption(optionData);
                    option->setPricingEngine(bsmhwEngine);
                    Real expected = option->NPV();

                    option->setPricingEngine(fdEngine);
                    Real calculated = option->NPV();
                    avgPriceDiff
                            += std::fabs(expected - calculated) / LENGTH(strikes);
                }

                if (controlVariate[i] && tolWithCV[l] < avgPriceDiff) {
                    FAIL_CHECK("Failed to reproduce BSM-Hull-White prices"
                                       << "\n   scheme       : " << scheme.name
                                       << "\n   model        : " << hestonModelData.name
                                       << "\n   CV           : on");
                }


                if (!controlVariate[i] && tolWithOutCV[l] < avgPriceDiff) {
                    FAIL_CHECK("Failed to reproduce BSM-Hull-White prices"
                                       << "\n   scheme       : " << scheme.name
                                       << "\n   model        : " << hestonModelData.name
                                       << "\n   CV           : off");
                }
            }
        }
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_SpatialDiscretizatinError", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing spatial convergence speed of Heston engine...");

    SavedSettings backup;

    Date today(27, December, 2004);
    Settings::instance().evaluationDate() = today;

    Real maturity = 1.0;
    Real strikes[] = {75, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 140, 150};
    Size listOfTimeStepsPerYear[] = {40};

    const Real tol[] = {0.02, 0.02, 0.02, 0.02, 0.05};
    for (Size u = 0; u < LENGTH(listOfTimeStepsPerYear); ++u) {
        for (Size i = 0; i < LENGTH(schemes); ++i) {
            for (Size j = 0; j < LENGTH(hestonModels); ++j) {
                Real avgPriceDiff = 0;
                std::shared_ptr < HestonProcess > hestonProcess(
                        makeHestonProcess(hestonModels[j]));
                std::shared_ptr < HestonModel > hestonModel =
                        std::make_shared<HestonModel>(hestonProcess);

                std::shared_ptr < PricingEngine > analyticEngine =
                        std::make_shared<AnalyticHestonEngine>(hestonModel, 172);

                Size tSteps = Size(maturity * listOfTimeStepsPerYear[u]);

                std::shared_ptr < FdHestonVanillaEngine > fdEngine =
                        std::make_shared<FdHestonVanillaEngine>(
                                hestonModel, tSteps, 200, 40, 0,
                                schemes[i].schemeDesc);
                fdEngine->enableMultipleStrikesCaching(
                        std::vector<Real>(strikes, strikes + LENGTH(strikes)));

                for (Size k = 0; k < LENGTH(strikes); ++k) {
                    VanillaOptionData optionData
                            = {strikes[k], maturity, Option::Call};
                    std::shared_ptr < VanillaOption > option
                            = makeVanillaOption(optionData);
                    option->setPricingEngine(analyticEngine);
                    Real expected = option->NPV();

                    option->setPricingEngine(fdEngine);
                    Real calculated = option->NPV();

                    avgPriceDiff
                            += std::fabs(expected - calculated) / LENGTH(strikes);
                }

                if (avgPriceDiff > tol[i]) {
                    FAIL_CHECK("\nFailed to reproduce Heston prices" <<
                                                                     "\n scheme    : " << schemes[i].name <<
                                                                     "\n model     : " << hestonModels[j].name <<
                                                                     "\n error     : " << avgPriceDiff <<
                                                                     "\n tolerance : " << tol[i]);
                }
            }
        }
    }
}


namespace {
    class HestonHullWhiteCorrelationConstraint : public Constraint {
    private:
        class Impl : public Constraint::Impl {
        public:
            explicit Impl(Real equityShortRateCorr)
                    : equityShortRateCorr_(equityShortRateCorr) {}

            bool test(const Array &params) const {
                const Real rho = params[3];

                return (square(rho)
                        + square(equityShortRateCorr_) <= 1.0);
            }

        private:
            const Real equityShortRateCorr_;
        };

    public:
        explicit HestonHullWhiteCorrelationConstraint(
                Real equityShortRateCorr)
                : Constraint(std::make_shared<HestonHullWhiteCorrelationConstraint::Impl>(
                equityShortRateCorr)) {}
    };
}


TEST_CASE("HybridHestonHullWhiteProcess_HestonHullWhiteCalibration", "[HybridHestonHullWhiteProcess]") {
    INFO("Testing the Heston Hull-White calibration...");

    // Calibration of a hybrid Heston-Hull-White model using 
    // the finite difference HestonHullWhite pricing engine
    //
    // Input surface is based on a Heston-Hull-White model with
    // Hull-White: a = 0.00883, \sigma = 0.00631
    // Heston    : \nu = 0.12, \kappa = 2.0, 
    //             \theta = 0.09, \sigma = 0.5, \rho=-0.75
    // Equity Short rate correlation: -0.5

    SavedSettings backup;

    const DayCounter dc = Actual365Fixed();
    const Calendar calendar = TARGET();
    const Date today = Date(28, March, 2004);
    Settings::instance().evaluationDate() = today;

    const Handle<YieldTermStructure> rTS(flatRate(0.05, dc));

    // assuming, that the Hull-White process is already calibrated 
    // on a given set of pure interest rate calibration instruments.
    std::shared_ptr < HullWhiteProcess > hwProcess =
            std::make_shared<HullWhiteProcess>(rTS, 0.00883, 0.00631);
    std::shared_ptr < HullWhite > hullWhiteModel =
            std::make_shared<HullWhite>(rTS, hwProcess->a(), hwProcess->sigma());

    const Handle<YieldTermStructure> qTS(flatRate(0.02, dc));
    Handle<Quote> s0(std::make_shared<SimpleQuote>(100.0));

    // starting point of the pure Heston calibration
    const Real start_v0 = 0.2 * 0.2;
    const Real start_theta = start_v0;
    const Real start_kappa = 0.5;
    const Real start_sigma = 0.25;
    const Real start_rho = -0.5;

    const std::shared_ptr<HestonProcess> hestonProcess =
            std::make_shared<HestonProcess>(rTS, qTS, s0, start_v0, start_kappa,
                                            start_theta, start_sigma, start_rho);
    const std::shared_ptr<HestonModel> analyticHestonModel =
            std::make_shared<HestonModel>(hestonProcess);
    const std::shared_ptr<PricingEngine> analyticHestonEngine =
            std::make_shared<AnalyticHestonEngine>(analyticHestonModel, 164);
    const std::shared_ptr<HestonModel> fdmHestonModel =
            std::make_shared<HestonModel>(hestonProcess);

    const Real equityShortRateCorr = -0.5;

    const Real strikes[] = {50, 75, 90, 100, 110, 125, 150, 200};
    const Time maturities[] = {1 / 12., 3 / 12., 0.5, 1.0, 2.0, 3.0, 5.0, 7.5, 10};

    const Volatility vol[] = {
            0.482627, 0.407617, 0.366682, 0.340110, 0.314266, 0.280241, 0.252471, 0.325552,
            0.464811, 0.393336, 0.354664, 0.329758, 0.305668, 0.273563, 0.244024, 0.244886,
            0.441864, 0.375618, 0.340464, 0.318249, 0.297127, 0.268839, 0.237972, 0.225553,
            0.407506, 0.351125, 0.322571, 0.305173, 0.289034, 0.267361, 0.239315, 0.213761,
            0.366761, 0.326166, 0.306764, 0.295279, 0.284765, 0.270592, 0.250702, 0.222928,
            0.345671, 0.314748, 0.300259, 0.291744, 0.283971, 0.273475, 0.258503, 0.235683,
            0.324512, 0.303631, 0.293981, 0.288338, 0.283193, 0.276248, 0.266271, 0.250506,
            0.311278, 0.296340, 0.289481, 0.285482, 0.281840, 0.276924, 0.269856, 0.258609,
            0.303219, 0.291534, 0.286187, 0.283073, 0.280239, 0.276414, 0.270926, 0.262173
    };

    std::vector<std::shared_ptr<CalibrationHelper> > options;

    for (Size i = 0; i < LENGTH(maturities); ++i) {
        const Period maturity(static_cast<int>(maturities[i] * 12.0 + 0.5), Months);
        std::shared_ptr < Exercise > exercise =
                std::make_shared<EuropeanExercise>(today + maturity);

        for (Size j = 0; j < LENGTH(strikes); ++j) {
            std::shared_ptr < StrikedTypePayoff > payoff = std::make_shared<PlainVanillaPayoff>(
                    strikes[j] * rTS->discount(maturities[i]) >=
                    s0->value() * qTS->discount(maturities[i])
                    ? Option::Call
                    : Option::Put,
                    strikes[j]);
            RelinkableHandle<Quote> v(std::make_shared<SimpleQuote>(vol[i * LENGTH(strikes) + j]));
            options.emplace_back(std::make_shared<HestonModelHelper>(maturity, calendar, s0,
                                                                     strikes[j], v, rTS, qTS,
                                                                     CalibrationHelper::PriceError));
            const Real marketValue = options.back()->marketValue();

            // Improve the quality of the starting point 
            // for the full Heston-Hull-White calibration 
            std::shared_ptr < SimpleQuote > volQuote = std::make_shared<SimpleQuote>();
            std::shared_ptr < GeneralizedBlackScholesProcess > bsProcess =
                    detail::ImpliedVolatilityHelper::clone(
                            std::make_shared<GeneralizedBlackScholesProcess>(
                                    s0, qTS, rTS, Handle<BlackVolTermStructure>(
                                            flatVol(v->value(), dc))), volQuote);

            VanillaOption dummyOption(payoff, exercise);

            std::shared_ptr < PricingEngine > bshwEngine =
                    std::make_shared<AnalyticBSMHullWhiteEngine>(equityShortRateCorr,
                                                                 bsProcess, hullWhiteModel);

            Volatility vt = detail::ImpliedVolatilityHelper::calculate(
                    dummyOption, *bshwEngine, *volQuote,
                    marketValue, 1e-8, 100, 0.0001, 10);

            v.linkTo(std::make_shared<SimpleQuote>(vt));

            options.back()->setPricingEngine(
                    std::shared_ptr < PricingEngine > (analyticHestonEngine));
        }
    }

    HestonHullWhiteCorrelationConstraint corrConstraint(equityShortRateCorr);
    LevenbergMarquardt om(1e-6, 1e-8, 1e-8);
    analyticHestonModel->calibrate(options, om,
                                   EndCriteria(400, 40, 1.0e-8, 1.0e-4, 1.0e-8),
                                   corrConstraint);

    options.clear();
    fdmHestonModel->setParams(analyticHestonModel->params());

    for (Size i = 0; i < LENGTH(maturities); ++i) {
        const Size tGrid = static_cast<Size>(std::max(10.0, maturities[i] * 10.0));
        std::shared_ptr < FdHestonHullWhiteVanillaEngine > engine =
                std::make_shared<FdHestonHullWhiteVanillaEngine>(fdmHestonModel, hwProcess,
                                                                 equityShortRateCorr,
                                                                 tGrid, 61, 13, 9, 0, true);

        engine->enableMultipleStrikesCaching(
                std::vector<Real>(strikes, strikes + LENGTH(strikes)));

        const Period maturity(static_cast<int>(maturities[i] * 12.0 + 0.5), Months);

        for (Size j = 0; j < LENGTH(strikes); ++j) {
            // multiple strikes engine works best if the first option
            // per maturity has the average strike (because the first option
            // is priced first during the calibration and the first pricing
            // is used to calculate the prices for all strikes
            const Size js = (j + (LENGTH(strikes) - 1) / 2) % LENGTH(strikes);

            std::shared_ptr < StrikedTypePayoff > payoff =
                    std::make_shared<PlainVanillaPayoff>(Option::Call, strikes[js]);
            Handle<Quote> v(std::make_shared<SimpleQuote>(vol[i * LENGTH(strikes) + js]));
            options.emplace_back(std::make_shared<HestonModelHelper>(maturity, calendar, s0,
                                                                     strikes[js], v, rTS, qTS,
                                                                     CalibrationHelper::PriceError));

            options.back()->setPricingEngine(engine);
        }
    }

    LevenbergMarquardt vm(1e-6, 1e-2, 1e-2);
    fdmHestonModel->calibrate(options, vm,
                              EndCriteria(400, 40, 1.0e-8, 1.0e-4, 1.0e-8),
                              corrConstraint);

    const Real relTol = 0.01;
    const Real expected_v0 = 0.12;
    const Real expected_kappa = 2.0;
    const Real expected_theta = 0.09;
    const Real expected_sigma = 0.5;
    const Real expected_rho = -0.75;

    if (std::fabs(fdmHestonModel->v0() - expected_v0) / expected_v0 > relTol) {
        FAIL_CHECK("Failed to reproduce Heston-Hull-White model"
                           << "\n   v0 calculated: " << fdmHestonModel->v0()
                           << "\n   v0 expected  : " << expected_v0
                           << "\n   relatove tol : " << relTol);
    }
    if (std::fabs(fdmHestonModel->theta() - expected_theta) / expected_theta
        > relTol) {
        FAIL_CHECK("Failed to reproduce Heston-Hull-White model"
                           << "\n   theta calculated: " << fdmHestonModel->theta()
                           << "\n   theta expected  : " << expected_theta
                           << "\n   relatove tol    : " << relTol);
    }
    if (std::fabs(fdmHestonModel->kappa() - expected_kappa) / expected_kappa
        > relTol) {
        FAIL_CHECK("Failed to reproduce Heston-Hull-White model"
                           << "\n   kappa calculated: " << fdmHestonModel->kappa()
                           << "\n   kappa expected  : " << expected_kappa
                           << "\n   relatove tol    : " << relTol);
    }
    if (std::fabs(fdmHestonModel->sigma() - expected_sigma) / expected_sigma
        > relTol) {
        FAIL_CHECK("Failed to reproduce Heston-Hull-White model"
                           << "\n   sigma calculated: " << fdmHestonModel->sigma()
                           << "\n   sigma expected  : " << expected_sigma
                           << "\n   relatove tol    : " << relTol);
    }
    if (std::fabs(fdmHestonModel->rho() - expected_rho) / expected_rho > relTol) {
        FAIL_CHECK("Failed to reproduce Heston-Hull-White model"
                           << "\n   rho calculated: " << fdmHestonModel->rho()
                           << "\n   rho expected  : " << expected_rho
                           << "\n   relatove tol  : " << relTol);
    }
}

TEST_CASE("HybridHestonHullWhiteProcess_H1HWPricingEngine", "[HybridHestonHullWhiteProcess]") {

    SavedSettings backup;

    /*
     * Example taken from Lech Aleksander Grzelak,
     * Equity and Foreign Exchange Hybrid Models for Pricing Long-Maturity
     * Financial Derivatives,
     * http://repository.tudelft.nl/assets/uuid:a8e1a007-bd89-481a-aee3-0e22f15ade6b/PhDThesis_main.pdf
    */

    const Date today = Date(15, July, 2012);
    Settings::instance().evaluationDate() = today;
    const Date exerciseDate = Date(13, July, 2022);
    const DayCounter dc = Actual365Fixed();

    const std::shared_ptr<Exercise> exercise =
            std::make_shared<EuropeanExercise>(exerciseDate);

    const Handle<Quote> s0(std::make_shared<SimpleQuote>(100.0));

    const Real r = 0.02;
    const Real q = 0.00;
    const Real v0 = 0.05;
    const Real theta = 0.05;
    const Real kappa_v = 0.3;
    const Real sigma_v[] = {0.3, 0.6};
    const Real rho_sv = -0.30;
    const Real rho_sr = 0.6;
    const Real kappa_r = 0.01;
    const Real sigma_r = 0.01;

    const Handle<YieldTermStructure> rTS(flatRate(today, r, dc));
    const Handle<YieldTermStructure> qTS(flatRate(today, q, dc));

    const Handle<BlackVolTermStructure> flatVolTS(flatVol(today, 0.20, dc));
    const std::shared_ptr<GeneralizedBlackScholesProcess> bsProcess =
            std::make_shared<GeneralizedBlackScholesProcess>(s0, qTS, rTS, flatVolTS);

    const std::shared_ptr<HullWhiteProcess> hwProcess =
            std::make_shared<HullWhiteProcess>(rTS, kappa_r, sigma_r);
    const std::shared_ptr<HullWhite> hullWhiteModel =
            std::make_shared<HullWhite>(Handle<YieldTermStructure>(rTS), kappa_r, sigma_r);

    const Real tol = 0.0001;
    const Real strikes[] = {40, 80, 100, 120, 180};
    const Real expected[LENGTH(sigma_v)][LENGTH(strikes)]
            = {{0.267503, 0.235742, 0.228223, 0.223461, 0.217855},
               {0.263626, 0.211625, 0.199907, 0.193502, 0.190025}};

    for (Size j = 0; j < LENGTH(sigma_v); ++j) {
        const std::shared_ptr<HestonProcess> hestonProcess =
                std::make_shared<HestonProcess>(rTS, qTS, s0, v0, kappa_v, theta,
                                                sigma_v[j], rho_sv);
        const std::shared_ptr<HestonModel> hestonModel =
                std::make_shared<HestonModel>(hestonProcess);

        for (Size i = 0; i < LENGTH(strikes); ++i) {
            const std::shared_ptr<StrikedTypePayoff> payoff =
                    std::make_shared<PlainVanillaPayoff>(Option::Call, strikes[i]);

            VanillaOption option(payoff, exercise);

            const std::shared_ptr<PricingEngine> analyticH1HWEngine =
                    std::make_shared<AnalyticH1HWEngine>(hestonModel, hullWhiteModel,
                                                         rho_sr, 144);
            option.setPricingEngine(analyticH1HWEngine);
            const Real impliedH1HW
                    = option.impliedVolatility(option.NPV(), bsProcess);

            if (std::fabs(expected[j][i] - impliedH1HW) > tol) {
                FAIL_CHECK("Failed to reproduce H1HW implied volatility"
                                   << "\n   expected       : " << expected[j][i]
                                   << "\n   calculated     : " << impliedH1HW
                                   << "\n   tol            : " << tol
                                   << "\n   strike         : " << strikes[i]
                                   << "\n   sigma          : " << sigma_v[j]);
            }
        }
    }
}
